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ABSTRACT 

Observations reveal concentrations of molecular line emission on the sky, called "clumps," in dense, 
star- forming molecular clouds. These clumps are believed to be the eventual sites of star formation. We 
study the three-dimensional analogs of clumps using a set of self-consistent, time-dependent numerical 
models of molecular clouds. The models follow the decay of initially supersonic turbulence in an isother- 
mal, self-gravitating, magnetized fluid. We find the following. (1) Clumps are intrinsically triaxial. This 
explains the observed deficit of clumps with a projected axis ratio near unity, and the apparent prolate- 
ness of clumps. (2) Simulated clump axes are not strongly aligned with the mean magnetic field within 
clumps, nor with the large-scale mean fields. This is in agreement with observations. (3) The clump 
mass spectrum has a high-mass slope that is consistent with the Salpeter value. There is a low-mass 
break in the slope at ~ 0.5M Q , although this may depend on model parameters including numerical 
resolution. (4) The typical specific spin angular momentum of clumps is 4 x 10 22 cm 2 s _1 . This is larger 
than the median specific angular momentum of binary stars. Scaling arguments suggest that higher 
resolution simulations may soon be able to resolve the scales at which the angular momentum of binary 
stars is determined. 

Subject headings: ISM: clouds, ISM: Molecules, Magnetohydrodynamics: MHD, Methods: Numerical, 
Stars: Formation 



1. INTRODUCTION 

Star forming molecular clouds are observed to be highly 
inhomogeneous in molecular line emission and continuum 
absorption studies. They contain relatively empty voids 
as well as dense "clumps" and denser "cores" within these 
clumps. One expects that stars will form in the densest 
parts of molecular clouds, and clumps are indeed corre- 
lated with young stars (Williams et al. 1994; Phelps & 
Lada 1997). Clumps (or cores) therefore appear to be the 
immediate precursors of clusters, small groups, or individ- 
ual stars in the interstellar medium; they constitute the 
initial conditions for star formation. 

Theoretical models of diffuse precursors to stars were, 
until recently, confined to static or quasi-static equilib- 
rium models. Mouschovias (1976) constructed numerical 
models of static, axisymmetric, magnetically supported 
self-gravitating equilibria (see also: Barker & Mestel 1996 
Li & Shu 1996; Galli & Shu 1993; Baureis et al. 1989 
Tomisaka et al. 1988a,b; Mestel & Ray 1985). These mod- 
els do not include the effects of turbulence. Quasi-static 
models, in which the effects of turbulence are included 
as an effective pressure, fixed by an assumed equation of 
state, have been constructed by Lizano & Shu (1989) (see 
also Myers & Fuller 1992; McLaughlin & Pudritz 1996; 
Curry & McKee 2000). 

Both the static and quasi-static approach have weak- 
nesses. The static, magnetically supported models predict 
oblate cores, flattened perpendicular to the mean magnetic 
field (but see Fiege & Pudritz 2000). Observations yield 
shapes apparently more consistent with prolate clumps 
(Myers et al. 1991; Ryden 1996), and little correlation be- 



tween field direction and clump shape (Heilcs et al. 1993). 
They also do not include the effects of turbulent velocities, 
which are directly observed through their contribution to 
the linewidth of molecular line transitions even at small 
spatial scales (e.g. Goodman et al. 1998). 

The quasi-static models, on the other hand, account for 
turbulent velocities by invoking a steady turbulent pres- 
sure. It not immediately apparent that supersonic tur- 
bulence like that observed in molecular clouds would act 
to prevent collapse, since the turbulence itself leads di- 
rectly to compressions that might accelerate collapse. Re- 
cent numerical simulations show, furthermore, that su- 
personic magnetized turbulence dissipates in a dynamical 
time (Stone et al. 1998; Mac Low 1999). If these simu- 
lations correctly represent conditions in molecular clouds 
then even an approximate steady state cannot be achieved 
unless momentum is continually injected into the cloud. 

Three-dimensional (3D), dynamic numerical models of 
molecular clouds have recently become practical (e.g., 
Vazquez- Semadeni et al. 1998 and the references therein; 
Stone et al. 1998; Mac Low 1999; Padoan & Nordlund 
1999; Pichardo et al 2000; Ostriker et al. 2001). They do 
not impose a high degree of symmetry on the cloud and 
allow a largely self-consistent treatment of turbulence, al- 
though due to numerical diffusion they have a much lower 
effective Reynolds numbers than real interstellar clouds. 
In this paper we will study the density concentrations that 
arise self-consistently from self-gravitating turbulence in a 
set of numerical models. 

Before we proceed, it is worth discussing the limita- 
tions of the model and the approximations that we use. 
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Our models integrate the equations of isothermal, com- 
pressible, self-gravitating, ideal magnetohydrodynamics 
(MHD). The isothermal approximation is used because re- 
alistic, time-dependent heating and cooling that treats the 
transfer of line radiation out of the cloud is, at present, 
out of numerical reach (but see Juvela et al. 2001). The 
model is fully self-gravitating, however, we use periodic 
boundary conditions in our solution of the Poisson equa- 
tion. This distorts the gravitational field from what it 
would be in an isolated cloud by introducing additional 
tidal forces (although real clouds arc not isolated, tidal 
forces are smaller). The numerical model is "ideal" in the 
sense that we do not include any explicit physical diffu- 
sion coefficients. All dissipation occurs in shocks, where 
it is captured by an artificial viscosity, or via numerical 
diffusion at the grid scale. Thus we do not model phys- 
ical diffusion processes such as ion- neutral drift. This is 
almost certainly important in the dense regions of clouds, 
but an accurate treatment of the ionization, particularly 
UV ionization, that would motivate dynamical studies of 
non-ideal clouds is also numerically out of reach. 

This paper is organized as follows: In §2 we give our 
definition of a "clump." In §3 we describe our numerical 
methods and tests and give a qualitative and quantitative 
description of the evolution of the runs. The clump mass 
spectrum is described in §4. The distribution of clump 
shapes in described in §5, the orientation of the clumps 
with respect to the local and global magnetic field is de- 
scribed in §6, and a discussion on the clump spin angular 
momentum is given in §7. A summary and comparison 
with earlier work is given in §8. 

2. CLUMPS: DEFINITIONS 

Our first task is to define a "clump." Any definition is 
somewhat arbitrary, much as the designation "mountain" 
1 is arbitrary in topography. The definition should be: 
(1) physically motivated; (2) related to definitions used by 
observers in molecular line studies of clouds; (3) simple 
to describe; (4) easy to apply. Ultimately we are moti- 
vated by the desire to characterize those dense regions in 
molecular clouds that turn into stars. 

Myers (1999) adopted the following definitions: a re- 
gion of size ^ 10 pc over which mean hydrogen number 
density tih 2 ~ 10 2 cm~ 3 is called a "complex" , a region 
of size ~ 1 pc over which <~ 10 3 cm~ 3 is termed a 
"clump" and a "core" refers to a region of size <~ 0.1 pc 
over which n# 2 ~ 10 4 cm -3 . A complex is well sampled 
by the CO J = 1 — > transition, a clump can be traced by 
the J = 1 — > line of 13 CO while a core is traced by the 
(J, K) = (1, 1) line of NH 3 . The choice of 10 4 cm" 3 as the 
fiducial density has several advantages, e.g., some easily 
observable molecular lines have effective critical density of 
this order. 

There are two frequently-applied algorithms for identify- 
ing clumps in molecular line data. Each begins with inten- 
sity data in two space plus one velocity coordinates (l—b—v 
space). The first algorithm, GAUSSCLUMPS, uses a fit- 
ting procedure that removes objects with a constrained 
gaussian shape from the data (e.g. Stutzki & Glisten 
1990); it is analogous to the CLEAN algorithm used in ra- 
dio interferometry. The second algorithm, CLUMPFIND, 

1 Merriam- Webster: "a landmass that projects conspicuously above its 



identifies intensity maxima in I — b — v space (Williams et 
al. 1994), then assigns each "pixel" (i.e., each channel in 
a spectrum) to be a member of the clump associated with 
the nearest intensity maximum. 

These algorithms all implicitly make the assumption 
that the velocity coordinate v is somehow related to the 
line of sight distance z, i.e. that widely separated material 
that happens to lie at the same velocity is not confused 
in projection. As we have previously shown using simu- 
lated data sets (Ostriker et al. 2001; see also Pichardo et al 
2000; Ballesteros-Paredes & Mac Low 2002), this assump- 
tion may frequently not be satisfied for moderate-density 
condensations, although for rarer high-density peaks it 
probably is. With a suitable source function and radia- 
tive transfer scheme, it is possible to generated simulated 
l — b — v data cubes for molecular transitions with varying 
critical densities, and thus to ascertain exactly when v is 
expected to be a valid surrogate for z. We defer consider- 
ation of this interesting and important problem to a later 
publication. 

For the present work, we take advantage of our direct 
access to all three spatial coordinates in simulated data to 
study clumps in the unprojected density field. We asso- 
ciate clumps with spatial maxima in the density in much 
the same way that CLUMPFIND associates clumps with 
maxima in the intensity. Our algorithm operates as fol- 
lows. Find all density maxima rj. Then associate every 
zone r that lies above a threshold density p t with the near- 
est r,. Each rj, together with the surrounding material 
that lies above the threshold density, is then a clump. 

This algorithm is easy to implement and explain. It has 
one parameter p t . A potential problem is that the den- 
sity may contain small variation near the density maxima 
(are the Grand Tetons one mountain, or three?). If these 
variations occur on the grid scale, they may be purely nu- 
merical in origin. We can test the degree to which these 
variations split up what would otherwise be single clumps 
by smoothing the density and picking r, to be maxima 
of the smoothed density. This introduces one additional 
parameter, the smoothing length A. 

What is the relationship of the 3D clumps described here 
with those seen by observers, in projection? As mentioned 
above, preliminary work (Ostriker et al. 2001) indicates 
that projection effects can be important. In each section 
below we will apply the same clump analysis to both the 
column density and the density as a way of evaluating the 
importance of projection effects. 

3. NUMERICAL METHOD, TESTS, AND RUN EVOLUTION 

We evolve the equations of isothermal, self-gravitating, 
compressible, ideal MHD using a numerical method based 
on the ZEUS code (Stone & Norman 1992a,b). It is a 
time-explicit, operator-split finite difference algorithm on 
a staggered mesh. Density and internal energy are zone- 
centered, while velocity components live on zone faces. 
The magnetic portion of the code uses the "method of 
characteristics" to evolve the transverse components of the 
magnetic and velocity field in a manner that assures the 
accurate propagation of Alfven waves. The magnetic field 
is represented by fluxes through zone faces; these fluxes are 
evolved using "constrained transport," (Evans & Hawley 

surroundings and is higher than a hill." 
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which preserves the constraint V ■ B = to machine 
precision. 

The gravitational acceleration is obtained by differenc- 
ing the gravitational potential. The potential is calculated 
using a Fourier transform technique. The gravitational 
kernel is defined by 0k = -inGpk/ Q3,-(2 cos(Av Aa: ; ) — 
2)/ (Axi) 2 ). This kernel ensures that the natural finite 
difference form of the Poisson equation is satisfied. Notice 
that in a periodic domain the Poisson equation becomes 



V 2 c 



^G(p-p). 



(1) 



Evidently the periodic boundary conditions alter the po- 
tential from what it would be in an isolated cloud. The 
difference is most pronounced away from density maxima. 
In regions with p^> p the potential is less affected by the 
periodic boundary conditions. This can be verified by ex- 
plicit construction of sheetlike equilibria, which we have 
done. 

We use an isothermal equation of state p = c 2 p, c 2 — 
const. This models the effect of efficient heating and cool- 
ing, and is appropriate when the cooling time is much 
shorter than the dynamical time. 

Numerical convergence is a serious concern for the calcu- 
lations presented here, particularly for the case of clumps 
that are only a few zones in size. The self-gravitating 
isothermal MHD equations will, under a broad range of cir- 
cumstances, evolve toward density singularities (i.e. stars) 
that are unresolved. Many aspects of "protostellar" tur- 
bulence, however, may be insensitive to the internal struc- 
ture of small, high-density regions. Certainly this is the 
assumption on which our numerical investigation is un- 
dertaken. At a minimum, however, we must only measure 
quantities on well-resolved scales. 

3.1. Initial and Boundary Conditions 

This paper considers a set of three numerical models 
which differ in magnetic field strength. Several aspects of 
evolution and structure for the same three models have 
previously been investigated by Ostriker et al. (2001). 
Each model has a numerical resolution of 256 3 zones. 

The models are evolved in a periodic cubic domain of 
length L. The initial density and magnetic field are uni- 
form. The initial velocity field, 5v, is chosen to be a 
divergence-free Gaussian random field with power spec- 
trum v 2 nu fc~ 4 . This slope is compatible with Larson's 
law; see, e.g., Myers & Gammie (1999). The perturba- 
tions have E w = \ ] d 3 x(pSv 2 + SB 2 /An) = 100 pL 3 c 2 s , 
although initially SB — 0. No further perturbations are 
added over the course of the evolution; these are decay 
simulations. 

An important parameter of the model is the strength of 
selfgravity. At fixed physical density and sound speed, this 
strength is directly related to the physical size of the model 
L. A suitable nondimcnsional parametrization is in terms 
of the "Jeans number" nj = L/Lj, where L 2 = ixc 2 s l(Gp). 
All models discussed here have nj = 3. 

A second parameter is the strength of the mean mag- 
netic field, (B x ), where () indicates a volume average. 
(B x ) cannot change over time because of the periodic 



boundary conditions. A sense of the dynamical impor- 
tance of the mean field may be had by first comparing 
va, x = (B x ) /y/iirp to the sound speed c s . Our models 
have va,x/c s = 1 ("weak" field), \/T0 ("moderate" field), 
and 10 ("strong" field). The strong and moderate field 
cases are most consistent with observations (e.g. Crutcher 
1999; but see Padoan & Nordlund 1999 for a different per- 
spective). 

One may also compare the mean field strength to 
that required to prevent gravitational collapse: (B x ) > 
2-KpL\fG (Tomisaka 1988a, b). In terms of our dimension- 
less parameters va, x /c s and nj, this becomes va,x/c s > 
nnj. Models which satisfy this condition are subcritical, 
i.e., will never experience complete gravitational collapse, 
while systems which fail to satisfy it are supercritical. 2 In 
the present work, the strong field model is subcritical and 
the other two models are supercritical. 




Fig. 1. — A color image of the surface density from snapshot B 
9. The scale bar on the left shows column density relative to the 
mean column density. 

3.2. Run Evolution 

Once the model evolution is started, the initially incom- 
pressive velocity field distorts the initially uniform mag- 
netic field. Since the magnetic field distortions are typi- 
cally large, in the sense that 5B 2 /(8ir) > pc 2 , they couple 
strongly to compressive motions, and density variations 
result. Density variations also result from exchange of 
power between the compressive and incompressive parts 
of the velocity field; this would occur even in the absence 
of a magnetic field. Because the initial power spectrum 
of velocities we have chosen has v\, the velocity varia- 
tion at a scale A, varying as A 1 / 2 , and the timescale for 
a density variation to form is \/v\, the density variations 
appear first at the smallest scales. Within a dynamical 
time, L/vl, density variations have appeared at all scales 
and the flow is fully turbulent. 



2 The mass-to-flux ratio in our models is always constant. Relaxing this constraint introduces free functions, rather than a free parameter, 
in the initial conditions. These degrees of freedom will be explored in future models. Notice that deciding whether a cloud is sub- or super- 
critical requires knowledge of the field strength and geometry. 
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Each run begins with a kinetic energy of 100 ph 3 c\, 
which is equivalent to an r.m.s. Mach number (M = 
<t v /c s ) of \/200. This kinetic energy is rapidly reduced by 
work done in bending the mean field and in compressing 
the isothermal fluid. It is also dissipated in shocks. The 
total "wave energy" (difference between the energy and 
that of the initial state in the absence of perturbations) 
dissipates on a timescale of order L/vl, or 0.07L/c s . Fig- 
ure 1 shows a color-coded image of the surface density in 
the middle of one run (B); here black is low density and 
red is high density. 

In this paper we study a set of three identical experi- 
ments that differ only in the strength of the mean field. 
These are designated runs A (strong field, /3 = <? s /v\ x = 

0. 01), B (moderate field, = 0.1), and C (weak field, 
(3=1). Model A is subcritical, and models B and C 
are supercritical. Runs B and C undergo gravitational 
collapse in the sense that at some locations the density 
becomes very large. In run A the peak density rises to 
about 100/5 after the simulation begins and then drops to 
only 40p at the end of simulation. In contrast, the peak 
density rises almost monotonically in runs B and C, and 
at final stages, reaches 7.9 x 10 4 p (B) and 4.1 x 10 4 p (C). 
Note that this density violates the Truclove criterion 3 
(Truelove et al. 1997). At the end of B, the fraction of the 
mass lying at densities greater than (1, 10, 10 2 , 10 3 , 10 4 )p 
is (75,20,4,2,1)%. At the end of C, the fraction of the 
mass lying at densities greater than (1, 10, 10 2 , 10 3 , 10 4 )p 
is (70, 13,3, 1.5,0.9)%. Thus while high density, collapsed 
regions do develop in both supercritical runs, only a small 
fraction of the mass is in the highest density regions. The 
fraction of mass we find in collapsing regions is consistent 
with observations of the efficiency of star formation in gi- 
ant molecular clouds as a whole (e.g., the 2% found by 
Myers et al. 1986). 

From each of the simulations we have picked a set of 
"snapshots" to study intensively. They are labeled ac- 
cording to the run and evolution time (see Table 1). Run 
B, the moderate field ((3 = 0.1) run, is our "fiducial" run, 

1. e. the one we believe is most similar to a molecular cloud. 
Snapshot B9 is our fiducial snapshot, since it is evolved 
enough that the velocity and magnetic field have had a 
chance to relax, but not so evolved that all turbulence has 
decayed away. 

3.3. Parameters 

We have to assign values to both the threshold density 
pt and the smoothing length A (given in units of the grid 
spacing) to obtain a population of "clumps" . In addition 
to these two parameters there are two scaling parameters 
that must be chosen to give a physical scale to the simula- 
tion (there are two rather than three because gravity sets 
a lengthscale in the model) . We fix the mean number den- 
sity uh 2 = 100 cm~ 3 and the temperature T — 10 K, corre- 
sponding to typical values found in molecular clouds. With 
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these scaling parameters p = nu 2 p = 
for p = 2.4 to p , c s = \JkT I p - 
nj c s \pK~l (Gp) = 6.45 pc (since nj 
mass M tot = pL 3 = 1532 M . 

Once the scaling parameters have been chosen there is 
a natural choice for the threshold density p t - the density 
corresponding to the critical density for CO (J = 1 — ► 0) 
line, ~ 10 3 cm~ 3 . 4 Since n H . 2 = 100 cm~ 3 , p t = 10 p. 

There is no equally simple way to set the smoothing 
parameter A. The purpose of smoothing the dataset be- 
fore picking out the density maxima is to avoid separa- 
tion of closely spaced density peaks into separate clumps; 
one does not want to attach too much significance to ob- 
jects separated by only one grid zone. This suggests fixing 
A ~ 1. The number of clumps decreases as smoothing 
is increased, however, so choosing A too large would re- 
sult in the loss of most clumps and most of the informa- 
tion contained in the dataset. We have experimented with 
A = 1, 1.5, and 2. All these give similar results. All results 
reported here are for A = 1.5 unless otherwise noted. To 
avoid making an arbitrary choice of this kind one may al- 
ternatively consider an ensemble of clump sets defined by 
allowing A to vary from 1 to N/2, where N is the number 
of zones. This is operationally similar to the technique 
employed in Ostriker et al. (2001); Ostriker (2002) (and 
also in the analysis of large-scale cosmological structure) . 




Fig. 2. — A contour plot of the surface density in snapshot B9 
from our fiducial run. 

Are the clumps well resolved? To decide this issue, we 
define an effective radius r 2 ^ = J pr 2 /M cl , where the in- 



3 This criterion checks if the size of the Jeans length is of order a zone size or smaller. The instantaneous Jeans length is X 2 T = nc 2 J(Gp). smcc 

1 C 2 iV 2 7T 

c s is constant in space and time, we can rewrite the Truelove criterion as p < -^j q l '2 i here N is the number of resolution elements along 
each axis, L is the size of the box, and T , the "Truelove number" must be larger than about 4. In simulations units this becomes p < 455p. 
Therefore we are not resolving the internal dynamics of the clump with the highest peak density. 

4 In a real cloud, of course, radiative trapping permits regions of density lower than the critical value to contribute significantly to the observed 
emission, so the 12 CO-emitting regions in a real cloud might be somewhat larger than our defined clumps; clumps observed in optically-thinner 
lines like 13 CO or C ls O clumps may be closer in size to the condensations defined by our density threshold. 
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Table 1 
Runs and Run Snapshots 
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tegral is taken over the clump volume, r is measured from 
the clump center of mass, and M c i is the clump mass. The 
distribution of r e // depends on A. For A = 1.5 the clump 
size distribution in snapshot B9 (our fiducial snapshot) 
is such that 97% of the clumps have an effective radius 
r e ff > 3 grid zones, 78% have r e // > 5 grid zones, and 
15% have r e ff > 10 grid zones. Similar results are ob- 
tained for other snapshots. Thus most clumps are (at this 
stage of the evolution) well resolved. 




0.2 0.4 0.6 0.8 1 



x/L 

Fig. 3. — The location of 3D clumps (for definition see text) in 
snapshot B9, projected into the same plane as Figure 2. The center 
of each circle is at the clump center of mass, and its radius is equal 
to the clump effective radius. 



Are the 3D clumps identified with this set of parame- 
ters at all related to 2D clumps (i.e. the clumps found 
by applying our clump finding algorithm to column den- 
sity maps)? Figure 2 shows the column density contours 
for a line of sight aligned with one of the simulation box 
axes, and Figure 3 shows the positions of 3D clumps pro- 
jected onto the same plane. The 3D clumps tend to be 
located in regions of large column density, but the mass 
and volume fractions of 3D clumps and 2D clumps differ 
sharply. The 3D clumps comprise 0.142 M tot and 0.008 L 3 ; 
in projection the 3D clumps cover 39% of the cloud. The 
2D clumps comprise 0.325 M to t and 0.526 L 2 . The column 
density is thus a good tracer of the dense portions of molec- 
ular clouds, but much of the contribution to each column 
may come from regions with density below threshold or 
even from regions spatially separated from 3D clumps (see 
Ostriker et al. 2001). 

4. CLUMP MASS SPECTRUM 

The simplest property of simulated clumps that we 
can study is their mass spectrum. We characterize the 
mass spectra by four numbers: the maximum clump mass 
M max , the minimum clump mass M min , the mass corre- 
sponding to the peak of the spectrum M pea k (defined as 
the maximum of dN / dlnM for data binned in In M) and 
the slope of the high-mass wing, x. The slope x is de- 
termined by binning the clumps in InM and performing 
a linear fit to dN/d\nM oc Af - ^ -1 ), weighted according 
to the square root of the number of clumps in a bin in 
the range M pea k < M < M max . Table 2 gives a sum- 
mary of the mass spectrum characteristics for each of the 
snapshots. 

In measuring x our intent is simply to characterize the 
shape of the mass spectrum in a way that is similar to 
that used by observers, not to suggest that the underlying 
spectrum is a power law. Other forms of the mass spec- 
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trum may well provide a better fit. Also notice that the 
measurement of x is noisy because of the small number 
of clumps in the high mass wing. It can change by ±0.3 
depending on the choice of bin size and location. 

The minimum clump mass M m i n may be influenced by 
numerical resolution. We do not have a resolution study 
in hand and so we cannot quantitatively evaluate this. 
The minimum clump mass permitted by our clump find- 
ing algorithm is the mass in a single zone that exceeds the 
threshold density. This is (10/256 3 )M tot = 9.1 x 1CT 4 M Q , 
which is much smaller than our typical M m i n . This is re- 
assuring, but because the smallest clumps typically consist 
of a small number of zones the reader is cautioned against 
placing too much weight on our measurements of M m i n . 

Consider the mass spectrum of the fiducial snapshot 
B9, which is in many respects typical (Figure 4, solid 
line). A total of 300 clumps were identified, account- 
ing for 14.2% of the total mass and 0.871% of the vol- 
ume. The slope of the spectrum is x — 2.35 (it is an 
accident that this is so close to the Salpeter slope for 
the stellar initial mass function [IMF], since different val- 
ues of the smoothing parameter A and different binning 
schemes give values of x that vary by as much as ±0.3). 
The maximum mass M max = 5 x 10~ 3 M to t, the mini- 
mum mass M min = 6 x 10~ 5 M to t, and the peak mass 
Mp ea k = 3 x 10~ 4 M to t- In physical units the peak mass 
is close to 0.5 M Q , and the minimum and maximum to 
0.1 and 8 M , respectively. Clearly this is an interest- 
ing range of masses if clumps like these are the immediate 
precursors of star formation, but it is important to remem- 
ber the numerical limitations of our study. In particular 
we do not know how these values might vary with initial 
conditions, Jeans number n.j, Mach number, or numerical 
resolution. 

We expect that some clumps are transients, the result 
of collisions between oppositely moving streams of gas in 
which gravity plays no role, while others are long lived and 
self-gravitating. To measure the importance of gravity for 
each clump we define 



a 



2E k + 3E P + E b ' 



(2) 



where E g = (1/2) J{p - p) <pd 3 r 5 , E k = (1/2) J p{v 



„2 j3 



n is the clump center-of-mass velocity; E p = 
J pc 2 d 3 r is related to the internal energy and Eb = 
J B 2 /{Sit) d 3 r is the magnetic energy inside the clump; 
in each case the integral is taken over the clump volume. 
We somewhat arbitrarily label clumps as self-gravitating 
if a > 0.5 and nonself-gravitating otherwise. Notice that 
our a is slightly different from the usual "virial param- 
eter", which is proportional to the ratio of gravitational 
potential energy to kinetic energy. We will analyze the 
kinematics of clumps by the virial theorem in greater detail 
in a later paper (see, e.g., Zweibel 1990; McKee & Zweibel 
1992; Bertoldi & McKee 1992; Ballesteros-Paredes et al. 
1999). 

We show in Figure 4 the mass spectra for self-gravitating 
(dotted line) and nonself-gravitating (dashed line) clumps 
in snapshot B9. At this stage only 4% of the clumps 
(about 14% of total clump mass) are self-gravitating, and 
all self-gravitating clumps are in the high-mass wing of 

5 This is the gravitational contribution to the total, conserved energy 



the spectrum. The majority of the clumps are nonself- 
gravitating, and thus x for all clumps and for the nonself- 
gravitating subset are about the same. 
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Fig. 4. — Mass spectra of self- and nonself-gravitating clumps and 
all clumps from snapshot B9. The mass unit is c 3 a G -1,5 p~ ' 5 ~ 
1O.2M . 

How does the clump mass spectrum evolve (see, e.g., 
Williams & Blitz 1998)? In Figure 5 we show the mass 
spectra from the fiducial, earliest and last snapshots in 
run B. Several features stand out: (1) M max increases 
with time. From the first to last snapshot M max has in- 
creased about an order of magnitude. On the other hand, 
M m i n remains about the same. (2) M pea k remains about 
the same. (3) The number of clumps at peak (per bin) 
and the total clump number decrease with time. This is 
a result of clump agglomeration. Notice that, although in 
runs B and C the total clump number decreases with time, 
the total mass contained in clumps increases. Evidently 
clumps also accrete mass from lower-density material. (4) 
The slope of the high-mass wing becomes "shallower" , i.e., 
x decreases: x(B3) = 3.15, x(B9) = 2.35, ir(B19) = 1.98. 
(5) The number of self-gravitating clumps increases, es- 
pecially in the high-mass wing: in the fiducial snapshot 
there are only 13 self-gravitating clumps, but in the last 
snapshot the entire high-mass wing is comprised of such 
clumps (extending down to M pea k ~ 0.6 M Q ). 

Docs the field strength influence the mass spectrum? 
Figure 6 shows mass spectra from snapshots A19, B9, and 
C9. We have chosen to compare these snapshots because 
all have similar sonic Mach numbers M. ~ 5. For the two 
supercritical runs (B and C) the mass spectra are similar. 
The subcritical run (A), however, is different: there are 
essentially no self-gravitating clumps in any of the snap- 
shots and the mass spectrum changes in time in such a 
way that the number of clumps decreases as the evolution 
progresses. 

under periodic boundary conditions. 
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Table 2 
Statistics of Run Snapshots 
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"Mass fraction in clumps. 

'Mass fraction in self-gravitating clumps. 

c The uncertainty in the mass spectrum slope x is ~ ±0.3 

ing clumps as "2D clumps". We find that (1) masses 
of 2D clumps are larger than the 3D clumps (by about 
an order of magnitude) and therefore, M max and M p are 
both greater than those in 3D spectra. (2) The slopes 
are generally slightly steeper than the 3D ones and have a 
weak tendency to become shallower in time. For example, 
x(B3,2L>) = 3.26, x(B9,2D) = 2.49, a;(B19,2L>) = 2.50. 
The 2D and 3D mass spectra for snapshot B9 are com- 
pared in Figure 7. 

5. CLUMP SHAPES 

Clump shapes may be an indicator of the dynamical 
processes that govern clump internal dynamics. For ex- 
ample, early models of magnetized, self-gravitating clouds 
(e.g. Mouschovias 1976) have oblate spheroidal isodensity 
contours that are flattened in a plane perpendicular to the 
mean magnetic field. Observers have looked for evidence 
of such oblate structures in molecular line maps of clouds. 

Myers et al. (1991)(see also David & Verschueren 1987; 
Fleck 1992; Curry & Stahler 2001) measured the projected 
axis ratios (r) of 16 dense cores in dark clouds and found 
that (r) = 0.5 — 0.6. They modeled ensembles of identical 
spheroidal cores, either all prolate or all oblate, assuming 
that their orientation in space is isotropic. Their results 
showed that the prolate ensemble can account for (r) with 
reasonable intrinsic axis ratio, while it requires extremely 
flattened oblate cores to produce same projected axis ra- 
tios. This suggests that prolate objects are preferred, con- 
trary to the prediction of static, magnetized models. 

Ryden (1996) obtained a tighter constraint on the 
shapes of clouds by studying the full distribution of ap- 
parent axis ratios f(r). For several samples of dense in- 
terstellar clouds, including Bok globules, dense cores, and 
clumps, she mapped the estimated distribution /(r) into 
a distribution of intrinsic axis ratios under the hypothesis 
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Fig. 5. — The evolution of the clump mass spectrum in the 
fiducial run B. 

Finally, do projection effects change the mass spec- 
trum? We have examined the mass spectrum of clumps 
in column density maps made from our simulations. We 
use precisely the same clump definition for the column 
density maps as in our 3D density fields, except that 
no smoothing is applied (A = 0, since the transforma- 
tion to column density is effectively a smoothing opera- 
tion) and the threshold column density is set to be 1.1 
times the mean column density. We refer to the result- 



8 



Gammie et al. 



that the clouds were oblate or prolate spheroids. She found 
that for most of her cloud samples the data were inconsis- 
tent, or at best marginally consistent, with the hypothesis 
that clouds are oblate spheroids because of a lack of almost 
circular (r ~ 1) objects. The samples were consistent with 
the hypothesis that the clouds were prolate spheroids. 
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Fig. 6. — The dependence of clump mass spectrum on initial mag- 
netic field strength in snapshots from each of the three runs. The 
snapshots are selected so that the turbulent sonic Mach number is 
nearly constant. 
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Fig. 7. — Comparison of clump mass spectra for clumps iden- 
tified in three dimensions (3D clumps) with those identified from 
column density maps in B9. The differences arise in part from the 
overlapping of 3D clumps in projection. 



We have performed a similar analysis of a subset of 
the catalog of clumps in the M17 SW cloud developed 
by Stutzki & Giisten (1990) using their GAUSSCLUMPS 
clump finder. We have restricted attention to those clumps 
with minor axis FWHM larger than 20" to minimize the 
effects of rounding by the 13" FWHM beam, and we 
have subtracted the beam size in quadrature from the axis 
lengths given in the paper. Notice that GAUSSCLUMPS 
uses a slightly different definition of clumps and princi- 
pal axis lengths from ours. Given this sample of clump 
principal axis lengths, it is possible to construct an esti- 
mate for the continuous distribution f(r) via the nonpara- 
metric kernel estimator used by Ryden (we use a kernel 
width h = 0.08, on a sample of 41 clumps). The result is 
shown in Figure 8 as a bold dashed line (Figure 9 shows 
the same result applied to "2D clumps"; see below). The 
light dashed lines show the uncertainties in f(r) estimated 
via a bootstrap resampling technique. The mean axis ratio 
(r) = 0.61. Recall that (r) > 0.5 for oblate spheroids. 

As for elliptical galaxies one can only infer the intrin- 
sic shape of a population of spheroidal clouds from a set 
of axis ratios. If the underlying population is triaxial, no 
such inference is possible, since a mapping from a distri- 
bution with one axis ratios to that with two axis ratios is 
degenerate. 
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Fig. 8. — The distribution of apparent axis ratios of 3D clumps 
(histogram) in our model and the distribution of apparent axis ra- 
tios in Stutzki & Glisten's (1990) catalog of clumps from M17SW 
derived using Ryden's (1996) method. 

5.1. Clump Shape Definition 

We characterize clump shapes using the eigenvalues of 
the moment of inertia tensor: 



d r pxi Xj. 



(3) 



Here x is measured relative to the clump center of mass. 
The eigenvalues are Ma 2 > Mb 2 > Mc 2 , where M is 
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the clump's mass and a, b, c are the principal axis lengths. 
The corresponding eigenvectors e a , e;,, e c are the principal 
axes. The intrinsic shape can be characterized via the axis 
ratios $ = b/a and 7 = c/a (note that (3 differs from the 
plasma j3). 

One advantage of using the moment of inertia tensor to 
measure clump shape is that it can be related to observa- 
tions, which are done in two dimensions on the sky. The 
projection of the moment of inertia tensor is 



P = 
1 v — 



(4) 



Here X\^ are Cartesian coordinates on the plane of the 
sky and £ is the surface density or anything propor- 
tional to the surface density, such as total flux in an opti- 
cally thin line that traces the density. P has eigenvalues 
Mp 2 > Mq 2 . The apparent shape is characterized by the 
apparent axis ratio r = q/p. This is what can be measured 
by an observer from column density maps of an isolated 
clump. 



3 - 



T 



1 1 1 | 1 1 1 1 1 1 1 1 T- 



observation 
uncertainties 
in observation 
model 

/ 



/ \ 



A 











\ 




\ \ 






\ 




\ 



0.2 



0.4 



0.6 



0.8 



Fig. 9. — The distribution of axis ratios for 2D clumps (his- 
togram) together with the M17SW clump distribution from Stutzki 
& Glisten's (1990), as in Figure 8. 

The apparent shape can be calculated from the eigenval- 
ues a, b, c of the moment of inertia tensor without actually 
projecting the clump pixel by pixel (this is computation- 
ally useful). Let 9 be the angle between e c and line-of- 
sight, and 4> be the angle between e a and x, which points 
along the intersection of the e a , e b plane with the plane of 
the sky. Then the apparent axis ratio is 

1 I \ V2 

Pxx + Pyy ~ <J Pxx + ^Pxy ~ ^PxxPyy + Pyy 

^ Pxx + Pyy + \J Pxx + ^P-ry ~ %PxxPyy + Pyy 

(5) 

where P xx — a 2 sin 2 9 + b 2 cos 2 9 sin 2 4> + c 2 cos 2 9 cos 2 (j>, 
Pxy = {b 2 — c 2 ) cos 9 cos 4> sin cf>, P yy — b 2 cos 2 4> + c 2 sin 2 cf>. 



Finally, some terminology: the term "projected clump" 
as used here refers to an object found by our clumpfinder in 
three spatial dimensions (a 3D clump) and then projected 
as if it were completely isolated onto the sky. The term 
"2D clump" refers to an object found by our clumpfinder 
in two spatial dimensions, after the entire density field 
has been projected on the sky. Were all 3D clumps well 
separated these procedures would yield identical results. 
However 3D clumps overlap in projection and therefore 
the populations differ. 
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Fig. 10. — The distribution of 3D clump axis ratios in B9. 
Clumps lying on the diagonal are prolate (b/a = c/a < 1); clumps 
lying along the right hand boundary of the plot are oblate (a/b = 1; 
c < 1). We have divided the diagram into three regions: prolate 
(within the triangle marked P); triaxial (within the triangle marked 
T); and oblate (within the triangle marked O). The precise bound- 
aries of these regions are arbitrary. 

5.2. 3D Axis Ratios 

Figure 10 shows the distribution of true axis ratios in 
B9. The self-gravitating clumps are denoted as solid dots, 
while nonself-gravitating clumps are crosses. The clumps 
are truly triaxial in the sense that they do not cluster 
around the prolate or oblate axes. As a convenience, we 
divide the f3 — 7 plane into three parts: a "prolate" part 
(which lies above the line connecting (/3, 7) = (1.0,1.0) 
and (0.33,0.0)), an "oblate" group (which lies below the 
line connecting (/3, 7) = (1.0,1.0) and (0.67,0.0)), and a 
"triaxial" group (everything else). Of the 300 clumps, 41% 
are prolate, 47% are triaxial, and only 12% are oblate. Al- 
though the relative clump numbers vary from snapshot to 
snapshot, in general the prolate and triaxial groups com- 
prise about 90% of the clumps. A purely spheroidal (oblate 
or prolate) model for the intrinsic clump shape distribu- 
tion is not viable. 

The prolate, oblate, and triaxial fraction for each snap- 
shot are listed in Table 3. The only trend worth remarking 
on is the growth in prolate clumps with time in run A. As 
we shall see below, these objects have a weak tendency 
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to line up so that the long axis is perpendicular to the 
magnetic field. 

The broad distribution of clump shapes in the — 7 
plane is likely telling us that clumps have complicated in- 
ternal dynamics- they are turbulent- and are not well de- 
scribed by highly symmetric equilibrium models. The self- 
gravitating clumps, however, are slightly "rounder" than 
nonself-gravitating ones; they tend to lie closer to the point 
(/3,-y) = (1.0,1.0) than their nonself-gravitating counter- 
parts (this is not a resolution effect, as self-gravitating 
clumps tend to be larger and hence better resolved than 
nonself-gravitating clumps). 

5.3. Apparent Axis Ratio 

The procedure for obtaining the apparent axis ratio of a 
3D clump was described in §5.1. If we apply this procedure 
to every clump and average over many viewing angles, we 
can obtain a mean, or expected distribution for f(r). Fig- 
ure 8 (solid line) shows the distribution of the apparent 
axis ratios obtained from B9, as well as the result from 
M17SW (Stutzki & Glisten 1990). 

Crudely speaking, the two distributions are consistent in 
the sense that both show a characteristic absence of nearly 
round objects. It is this absence that most strongly rules 
out a parent population of oblate spheroids (Ryden 1996). 
The f(r) changes with time in the sense that the mean 
apparent axis ratio increases. In fact, /(r) for snapshots 
Bll and B19 lies within the error bounds permitted by 
the observations. A similar increase in the mean apparent 
axis ratio with time is seen in run C (the clumps become 
rounder with time), but not run A. 

Our discussion up to now has concerned apparent axis 
ratios of individual clumps. Clumps exhibit considerable 
overlap ("confusion"), however, so the distribution of axis 
ratios of apparent (2D) clumps may be different from the 
distribution of apparent axis ratios of real (3D) clumps 
seen in projection. Figure 9 shows a comparison between 
the former (2D clumps obtained from a single line of sight 
projection) and the M17 observations. Broadly speaking, 
they are consistent with each other and the 2D clumps 
also show the absence of round clumps. 

6. CLUMP ORIENTATION 

Steady-state theoretical models of magnetically sup- 
ported clouds suggest that clouds should be flattened per- 
pendicular to the magnetic field (e.g., Mouschovias 1976; 
Tomisaka et al. 1988a, b; Baureis ct al. 1989; Lizano & Shu 
1989; Galli & Shu 1993; Fiedler & Mouschovias 1993; Li & 
Shu 1996) or, in the language of this paper, that the clump 
minor axis should be nearly parallel to the local field. Ob- 
servations by Goodman et al. (1990) (also Goodman et al. 
1992) suggest, however, that cloud minor axes are not cor- 
related with the local field direction. Recent observations 
showing alignment (Matthews & Wilson 2000) arc for a 
small number of clumps in individual clouds and probably 
do not indicate a statistically significant alignment 6 . Here 
we examine correlations between the clump principal axes 
and the density-weighted mean field within clumps, as well 
as the mean global fields, (B). 

6 Assessing the statistical significance of alignments is complicated by the fact that both magnetic fields and clump orientations are spatially 
correlated. 



Fig. 11. — The distribution of angles between the shortest body 
axis and the density weighted mean magnetic field in the clump in 
B9. The dashed line shows the distribution expected if both vectors 
are chosen randomly. 



Figure 11 shows the distribution of the dot product of 
the clump short axis unit vector (e c ) and the unit vector 
associated with the density weighted mean magnetic field 
in the clump (hi) in snapshot B9. The dashed line shows 
the distribution expected if both vectors are chosen at ran- 
dom. Evidently there is no strong correlation between the 
mean local field and the clump minor axis. 

Figure 12 shows a similar plot for 3D clumps in projec- 
tion, also for B9. The distribution is obtained as follows: 
for a fixed viewing angle, we use the same procedure de- 
scribed in §5 to obtain the principal axes of the projected 
clump, as well as the density-weighted magnetic field vec- 
tor within each clump. Averaging over 10 4 viewing angles 
gives the distribution. Again there is no strong correlation 
between the local mean field and the clump minor axis. 

A more complete description of alignment statistics is 
given in Table 3. There we list, for each of the snapshots, 
the fractional alignment of clump axes and the global and 
local mean fields. The quantity C.A., for example, is the 
fraction of all clumps with |e n • b;| > 0.86, i.e. the fraction 
of all clumps with major axis closer than 30° to the local 
mean field. Similarly, CB. is the fraction of all clumps 
with intermediate axis closer than 30° to the local mean 
field, etc. 

Table 3 also describes the alignment of clump axes with 
the global volume mean magnetic field with unit vector 
h g . The quantity Q.A. is the fraction of all clumps with 
|e a • b s | > 0.86. The only significant global alignment is 
in the late stages of the strong field run, when the short 
axis has a weak tendency to line up with the global mean 
field. 
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Table 3 

Shape and Magnetic Field Alignment in Clump Snapshots 



local alignment global alignment proj. align. 
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0.10 


0.12 


0.09 


0.18 


0.36 


0.31 


Cll 


0.42 


0.45 


0.13 


0.14 


0.14 


0.12 


0.10 


0.13 


0.14 


0.36 


0.31 


C15 


0.53 


0.37 


0.10 


0.13 


0.14 


0.10 


0.10 


0.14 


0.18 


0.34 


0.31 



"Fraction of prolate (V), triaxial (T) and oblate (0) clumps. 

^Fraction of clumps whose density-weighted mean magnetic field is within 30 deg of the clump major axis (C.A-), within 30 deg of the intermediate 
axis (C.B.), or within 30 deg of the major axis (C.C.). 
c Fraction of clumps for which the global mean magnetic field is within 30 deg of the clump minor axis (Q.A..), etc.. 

d Fraction of clumps for which the projected density weighted mean magnetic field is within 30 deg of the projected minor axis (P.A.), etc.. 

iV is the number of clumps (and similarly for C.B., C.C., 
etc.). This naive estimate for the dispersion in alignments 
is not accurate because clump orientations are probably 
spatially correlated (i.e. not independent). For exam- 
ple consider run C, whose weak initial field makes any 
global alignments physically implausible. Nevertheless, C 
exhibits global alignments that are "statistically signifi- 
cant" according to our naive estimate. 

The only significant and physically plausible alignments 
occur in run A, where we find that the minor axis tends 
to be weakly aligned with the field. The longest axis also 
tends to be anti-aligned with the field. This is similar to 
what might be expected from an equilibrium model, ex- 
cept that the clump shapes tend to be dominantly prolate 
or triaxial rather than oblate. 

Table 3 also describes the alignment of projected mean 
field vectors with the principal axes of the projected mo- 
ment of inertia tensor. These alignments are given in the 
final two columns (V.A. and V.B., which are the fraction 
of clumps with projected field less than 30° from the pro- 
jected clump major and minor axes, respectively). For 
randomly oriented clumps we expect V.A. = V.B. = 0.33. 
The projected alignments are sensitive to viewing angle, so 
the values shown have been averaged over 10 uniformly 
distributed viewing angles. The only plausible alignments 
are in run A; other snapshots show at best marginally 
significant alignments. To sum up, our numerical mod- 
els suggest that projected alignments are weak and not a 
useful diagnostic of field strength. 

Clump principal axis alignments, like the mass spec- 
trum, are sensitive to the smoothing parameter A. For 
example, if A = we see an implausible set of alignments 
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Fig. 12. — The distribution of angles between the projected 
shortest body axis and the projected density weighted mean mag- 
netic field in B9. The dashed line shows the distribution expected 
if both vectors are chosen randomly. 



Notice that if clump magnetic fields are randomly ori- 
ented with respect to the clump principal axes, and if 
the orientation of each clump is independent, we expect 
C.A. = 0.14±0.52/\//V at the 87% confidence level, where 



12 



Gammie et al. 



in run C. Presumably this is because we are studying align- 
ments near the grid scale. This is one of the motivations 
for introducing a finite A. 

7. ROTATION 

Observations suggest that rotation is not important in 
supporting clumps against gravitational collapse since ro- 
tational energy is small compared to the gravitational po- 
tential energy (e.g. Goodman et al. 1993). However, once 
collapse takes place, angular momentum may play an im- 
portant role in determining the characteristics of the sys- 
tem to be formed, e.g., the size of circumstellar disks, or 
the separation of the stars in a multiple star system. 

Clumps or cores are observed to have greater specific an- 
gular momentum j than binary or multiple star systems. 
Consider a binary with masses of 1 M© for both members 
and a period of 100 years. The specific angular momen- 
tum of this system is jbin = 1.16 x 10 20 cm 2 s -1 . This is 
typical for young binary T Tauri stars (e.g. Simon et al. 
1995). For single protostars and the envelope, Ohashi et 
al. (1997) obtained j ~ 10 20 cm 2 s -1 for two IRAS sources 
(for circumstellar disks, see, e.g. Beckwith et al. 1990). 
As for dense cores, the results given by Goodman et al. 
(1993) range from 6 x 10 20 to 4 x 10 22 cm 2 s" 1 . 

10' 10 2 j/Jwn io 3 io 4 
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Fig. 13. — Distribution of specific clump spin angular momen- 
tum in B9. The normalization is to c 2 G~ 0,5 p~°' 5 = jtot/(nj it ' 5 ), 
where j tot = 3.8 X 10 23 cm 2 s -1 = 3280 j bin and j bin a 10 20 cm 2 s -1 
is a reference angular momentum for a pair of 1 Mq stars in circular 
orbits at 100 AU. 

We have examined the specific angular momenta of the 
simulated clumps, and find that the values are greater than 
jbi n . In Figure 13 we show the distribution of specific an- 
gular momenta of clumps in B9. Clump specific angular 
momenta span about 2 — 3 decades, e.g., from 10 jbin to 
~6x 10 3 jbin, the peak is at 2.5 x 10 2 jbin', the range of 
the distribution decreases with time in B and C snapshots 
(we cannot say much about run A because there are not 
many clumps in later snapshots). There is no clear trend 



in clump specific angular momentum with mass. 

Our distribution of specific angular momenta may be 
compared to those found by Burkert & Bodenheimer 
(2000) (hereafter BB), who studied clump models in which 
the internal velocities are drawn from a Gaussian ran- 
dom field at a spatial scale ~ O.lpc; their power spectrum 
(v 2 oc fc -4 ) was normalized to match the typical velocity 
dispersion in the Goodman et al. (1998) sample of cores. 
Their distribution peaks at ~ 2 x 10 21 cm 2 sec -1 , while 
ours peaks at 4 x 10 22 cm 2 sec -1 . This difference is at- 
tributable to the following factors. First, our sample con- 
tains both nonsclf-gravitating and self-gravitating clumps, 
while BB's clumps are scaled to match self-gravitating 
clumps. Second, the self-gravitating clumps in our sample 
are more massive (by a factor of about 4) than BB's using 
our standard scaling. Third, the self-gravitating clumps 
have lower mean specific angular momentum (by about 
40%) than the nonsclf-gravitating clumps. Fourth, while 
BB chose a Gaussian clump density profile, our clumps 
have substantial envelopes that can make a significant con- 
tribution to clump angular momentum. We conclude that 
the difference is not significant. 

8. SUMMARY 

The main results of this paper are: 

1. Clumps are triaxial. Clumps that are "prolate" or 
"triaxial" (in the sense defined in §5) make up 90% of 
the clump population. Clumps that are "oblate" make up 
only 10% of the clump population. The distribution of 
clump shapes implies a distribution of apparent axis ra- 
tios that is consistent with observations. We do not see 
any signs of the sort of oblate equilibria that have been 
the focus of earlier studies of magnetized clouds emerging 
in our simulations. But we have not considered the most 
general possible set of initial states (we assume the mass- 
to-flux ratio is uniform), nor have we included the effects 
of ambipolar diffusion. 

The possibility of triaxial cloud cores has also been ex- 
plored by Basu (2000); Jones et al. (2001); Jones & Basu 
(2002). Basu (2000), however, suggested that the short 
axis should be aligned with the magnetic field. This is 
not seen in our numerical models. The absence of mag- 
netically supported equilibria in self-gravitating, turbulent 
MHD models has also been noted by several other work- 
ers in the field, perhaps most forcefully by Padoan et al. 
(2000). The emergence of such equilibria in supersonic, 
self-gravitating turbulence cannot yet be ruled out, how- 
ever, because the effects of ambipolar diffusion, varying 
mass to flux ratios, and varying initial conditions have not 
yet been studied. 

2. Clump principal axes are not strongly aligned with 
the magnetic field, either locally or globally. We would 
not expect to see such alignments in observations, nor are 
they currently observed. The effect of magnetic fields on 
cloud structure has been considered earlier by Ballesteros- 
Paredes & Mac Low (2002). They found that only strong 
fields produce an observable effect on the density structure 
based on a qualitative analysis. A similar result was seen 
at late stages of self-gravitating evolution in the strong- 
field models of Ostriker et al. (1999). 

3. The simulated clump mass spectrum has a character- 
istic mass of ~ 0.5 M© (for an assumed mean cloud density 
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and temperature). While this result is suggestive, it de- 
pends on some of the choices made in setting up our simu- 
lations, including numerical resolution. It is possible that 
the characteristic mass scales with the total cloud mass or 
with the sonic Mach number or density of the initial con- 
ditions. These scalings have not yet been explored. The 
high-mass wing of the clump mass spectrum has a slope 
that is crudely consistent with the Salpeter law, although 
other, non-power-law forms for the mass spectrum may 
well be consistent with our data. 

Other numerical work on the clump mass spectrum by 
Padoan ct al. (2000); Klessen (2001); Klessen & Burkert 
(2001); Padoan & Nordlund (2002); Ballesteros-Paredes & 
Mac Low (2002) is generally consistent with our results. 
These models include similar physics to ours, although 
some models lack selfgravity (Ballesteros-Paredes & Mac 
Low 2002) , while others lack magnetic fields (Klessen 
2001; Klessen & Burkert 2001). 

The simulations of Padoan et al. (2000) are the most 
similar to those analyzed here. Differences in Padoan et 
al.'s analytical methods limit the possibility of detailed 
comparisons, however: their definition of a self-gravitating 
clump does not include turbulent kinetic energy; the mass 
spectra presented are composites of two models with gas 
temperatures that differ by a factor of three; their defini- 
tion of a clump depends on contrast of local peaks rather 
than a density threshold. Nevertheless, the results of the 
two analyses for the position of the peak of the mass spec- 
trum, as well as for the slope of the high-mass end, are 
broadly consistent. 

Klessen (2001) and Klessen & Burkert (2001) use a com- 
pletely different numerical method (smoothed particle hy- 
drodynamics) and include selfgravity but not magnetic 
fields. In addition, they allow dense regions to turn into 
"sink" particles. This allows them to continue the evolu- 
tion much further in time than is possible in our calcula- 
tions. The trends we find of an increasingly flat high-end 
spectrum, and a larger proportion of mass in bound ob- 
jects, are clearly evident in their distributions at sequential 
evolutionary stages. Their clump mass spectra are again 
crudely consistent with ours; minor differences in mass 
scale and shape may arise due to differences in initial con- 
ditions or input parameters (number of Jeans masses). We 
cannot evaluate this without a direct attempt to reproduce 
their models. 

Recently, Padoan & Nordlund (2000) (see also Padoan 
& Nordlund 2002) have proposed that the clump mass 
spectrum, and consequently the mass function of collaps- 
ing cores, may be predicted solely from three macroscopic 
quantities averaged over the entire cloud: density, tem- 
perature, and velocity dispersion. Although this theory 
is attractive in many way, the stationary turbulent state 
adopted for the Padoan & Nordlund (2000) may not be a 
good approximation for real star forming regions. Clouds 
are transient structures, and a cloud may have developed 
its current mean density, temperature, and velocity dis- 
persion in many ways. Our analysis of purely decaying 
models shows that some aspects of the mass spectrum are 
invariant in time (e.g. M pea fc) even as the Mach number 
decreases, while other aspects of the mass spectrum (the 
fraction of self-gravitating clumps) increase secularly with 
time. Both of these results suggest that a cloud retains a 



memory of its dynamical history that may be important 
for establishing the stellar IMF. 

4. Clump specific angular momenta are an order of 
magnitude larger than in a "typical" binary system. It 
is interesting that our simulations give clump angular mo- 
menta that are so close to the angular momenta of binary 
systems. This suggests that we are close to resolving the 
lengthscale where the angular momentum of binaries is de- 
termined (since Av\ ~ A 1 / 2 in molecular clouds and in su- 
personic turbulence generally, the specific angular momen- 
tum j ~ Av\\ <~ A 3 / 2 , so A ~ j 2 / 3 )- Future higher reso- 
lution numerical models, perhaps with an adaptive mesh, 
will resolve this scale and may be capable of directly sim- 
ulating binary formation. 

Our results all depend on how the clumps are defined. 
Because fragmentation during gravitational collapse may 
decide the final stellar mass function, characterizing the 
hierarchy of clumpy structure and understanding its con- 
sequences may be very important. Mass functions are 
steeper when clumps-within-clumps are counted rather 
than discounted (see e.g. Ostriker 2002); whether gravity 
eventually selects the larger or smaller scale of a nested 
structure depends on the details of the collapse. More 
general tools (e.g. wavelet methods) for analyzing the den- 
sity field will be valuable in characterizing this hierarchi- 
cal structure. Future analyses should also focus on direct 
analogues of observable clumps (e.g. simulated clumps 
in position-position-velocity space for common molecular 
tracers), although self-consistent treatments of thermody- 
namics, radiative transfer, and molecular chemistry may 
ultimately be necessary for realistic source functions. 

This work was supported by NASA grants NAGW 5- 
9180 and NAG 59167. 
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